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We introduce and study a new directed sandpile model with threshold dynamics and stochas- 
tic toppling rules. We show that particle conservation law and the directed percolation-like local 
evolution of avalanches lead to the formation of a spatial structure in the steady state, with the 
density developing a power law tail away from the top. We determine the scaling exponents charac- 
terizing the avalanche distributions in terms of the critical exponents of directed percolation in all 
dimensions. 
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Many extended slowly driven dissipative systems in 
nature evolve into self-organized critical (SOC) steady 
states which show long-range spatial and temporal cor- 
relations. Since the pioneering work of Bak et al. [Q, 
sandpile models have served as paradigms of SOC sys- 
tems. A great deal of understanding of SOC has been 
achieved by numerically studying sandpile models with 
different evolution rules |||. For a special subclass of 
these models, i.e., the Abelian Sandpile Models (ASM), 
several analytical results are available |^ . Recent studies 
of models with stochastic toppling rules have shown that 
these models usually belong to a universality class differ- 
ent from deterministic automata; they have robust criti- 
cal states with respect to changes of a control parameter, 
and may also exhibit a dynamical phase transition be- 
tween qualitatively different steady states However, 
the spatial structures in the steady states of these models 
are much less investigated |^ . Due to long-range correla- 
tions in the critical states, influence of the boundary can 
be felt deep inside, and this can give rise to large-scale 
spatial structures. Indeed, emergent spatial structures 
are sine qua non for a SOC theory of fractals occurring 
in nature, e.g., mountain landscapes, river networks, and 
earthquake fault zones. 

In this Letter, we propose a new stochastic sandpile 
model which shows emergent spatial structures in the 
steady states. The model is a stochastic generalization 
of the directed ASM [|| and contains a probabilistic con- 
trol parameter p. In the case p ~1 the exponents are 
known exactly in all dimensions We show that for 
p 1, the model is in a new universality class and can 
be related to directed percolation (DP) ||^ problem with 
a nonuniform concentration profile. A steady state exists 
only for p > p* , where p* equals the critical threshold for 
directed site percolation. Above p* the system evolves 
towards a steady state which is arbitrarily close to the 
generalized DP critical line. We also show that in the 
critical state of our model a spatial structure results from 
an interplay of DP-like local evolution rules on the one 
side, and the dynamic conservation law on the other. We 



find a power-law density profile, which further enables us 
to determine the exact expressions for the scaling expo- 
nents of avalanches in terms of the DP critical exponents 
in all dimensions d. Our numerical simulations in two 
dimensions support these conclusions. 

For concreteness, we consider a square lattice of linear 
size L with sites (i,j) oriented so that the diagonal di- 
rection (1, 1) is vertically down. A non-negative integer 
variable, height h{i,j), is attached to each lattice site. 
Sand grains are added one at a time at a randomly cho- 
sen site on the top layer, increasing its height by one. A 
site becomes unstable if h{i,j) > 2 and relaxes as follows: 
With probability p, the local height decreases by two, and 
heights at each of its two downward neighbors increase 
by one. Otherwise, the heights remain unchanged. In ei- 
ther case, the site is considered as stable at the next time 
step. Only sites to which at least one particle was added 
at the preceding time step are checked for toppling. A 
discrete-time parallel update is applied to all such sites. 
We apply periodic boundary conditions in the horizon- 
tal direction. Two particles leave the system for each 
toppling occurring at the bottom layer. 

In the limit p ~ 1, the structure of the steady state 
is known exactly Only configurations with heights 
h = and h = 1 occur, and all such configurations are 
equally probable. The avalanche clusters are compact. 
For stochastic toppling, p ^ 1, the model has quite a 
different behavior. It is no longer Abelian, since adding 
two particles together to a site of height h > 2 does not 
have the same effect as adding them one by one at dif- 
ferent time steps. (In the latter case, toppling can occur 
twice with a finite probability. ) With a small probability 
heights could become arbitrarily large. The avalanche 
clusters have branches and holes of various sizes. An 
example is shown in Fig. 1. 

For small p the average influx of particles per at- 
tempted toppling at a site, which is > 1, exceeds the 
average outflux, 2p, thus there exists a value p*, such 
that a steady state is possible only for p > p* . We now 
argue that p* = Pc, where Pc is the critical threshold for 
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the directed site percolation problem 

Suppose that C is a stable configuration of the pile, 
and a particle added at the site (i, j) causes on the aver- 
age ni{C) topplings at depth £ below it. If for any stable 
configuration C , all heights are greater than or equal to 
corresponding heights in C then ne{C') > ne{C). Now, 
note that for all configurations C for which all sites have 
h > 1, the distribution of size of avalanches is exactly 
the same as in the directed sife-percolation problem on 
this lattice. Therefore, for all p < Pc, ne decreases ex- 
ponentially with £. Hence no topplings occur at large 
depths, and particles pile up in the upper layers, and 
thus there is no steady state. Conversely, for p > pc, the 
directed percolation avalanche clusters typically form a 
wedge, and ni increases with £ for large £ (see below). 
Then avalanches in configurations with all heights h > 1 
cause many topplings, and each layer after the avalanche 
on the average loses particles. Thus if the system ever 
reaches a state with large density, the number of parti- 
cles in the system will decrease until at sufficiently many 
sites heights become low enough so that the propagation 
of avalanches is affected, and it becomes critical, but not 
super- critical. Hence the system will have a steady state 
for all p > p* — Pc- On the square lattice numerical 
estimate for p^ § gives p* « 0.7054853(5). 

Our numerical simulations support this conclusion. In 
Fig. 2 the probability P{T) that an avalanche has dura- 
tions > T is plotted against T for different values of p 
and L ==200. (Notice that in directed models T = £). 
For large T it varies as P{T) ^ T^^'^K Exponential de- 
cay of the lower three curves indicates loss of SOC. In 
the steady states for p > p* the integrated cluster-size 
distribution behaves as D{s) ^ 5^^"^% with the following 
scaling properties 

D{s,L) = L^-^*V{sL-^^^), (1) 

and the scaling relation (r^ — 1) D\\ = — 1 is satisfied. 
In Fig. 3 the distribution D{s, L) vs. s and its finite-size- 
scaling plot are shown for p = p*. 

We now discuss the structure of the steady state for 
p* < p<l. Let p be the probability that a site chosen at 
random has a nonzero height in the steady state. Then 
the probability that this site will topple if a single par- 
ticle is added to it is Pi = pp, and the probability that 
it topples if two particles are added to it is P2 = p. The 
correlations of heights on the same layer are irrelevant 
and can be ignored Therefore, the evolution of an 
avalanche is the same as in a Domany-Kinzel (DK) cel- 
lular automaton model of directed site-bond percolation 
pof . This implies that, in order for the system to have 
critical correlations in the bulk, the set of points (Pi, P2) 
must lie on the critical line of the DK model. 

However, as we show below, the dynamic conservation 
law prevents the avalanche clusters in the steady state 
from being in the universality class of DP. Particle con- 
servation implies that in the steady state the average 



number of topplings at each layer equals 1/2. On the 
other hand, in the case of DP the expected number of 
growth sites at depth £ is known to vary as m with 
K = [{d — l)j^j_ — 2/3]/z^||, where /3, and are stan- 
dard DP critical exponents of the order parameter, and 
parallel and transverse correlation lengths, respectively. 
For DP in d<5 dimensions k>0, thus m increases with 
£, which is clearly not possible in the steady state. The 
way these conflicting requirements of particle conserva- 
tion and locally critical DP-like evolution are satisfied 
in our model is that the critical steady state develops a 
spatial structure. The density p, and hence Pi, are not 
uniform throughout the system, but vary from layer to 
layer |l|]. Let p{£) be the average density of sites with 
non-zero height in the £-th layer. By equating average 
infiux and outflux of particles at a site on ^-th layer, we 
find that p(£) = [l-(2p-l)/(^)]/[2p(l-/(£)], where /(£) 
is the number of topplings caused by simultaneous addi- 
tion of two particles at the site. The exactly calculated 
values of f{£) for the first few layers indicate that p{£) in- 
creases with £. As discussed above, for large £ the profile 
reaches the value p*{p) = P*{p)/p, where {Pi{p),p) is a 
point on the DP critical line in the (Pi,P2) parameter 
space. In Fig. 4, we plot the profile p against £ obtained 
by numerical simulations for p = p* and L — 200. The 
profile is well described by a power law: 

p{£)^p*{p)-A{p)£-'^ , (2) 

with p* =1 and A=G.39 for p = p*, and x ==0.578. 

We now show that the profile given by Eq. (||) changes 
the avalanche statistics in our model, and thus strongly 
affects the bulk transport. In the bulk, transport of par- 
ticles is locally described by the DK model with param- 
eters {Pi{£),P2). For large £ the system is close to the 
critical line and the local longitudinal correlation length 
^(£) varies as £,{£) - [Pj* - Pi (€)]"''" . In order for the 
transport to propagate further to a distance £, £,{£) must 
be proportional to £, i.e., 

m « t^/B . (3) 

This implies that the exponent x in Eq. (^ is exactly 
X — From the simulation data in log(l — p{£)) vs. 

log^ plot wc find the slope x = 0.575 ± 0.005 (see in- 
set to Fig. 4), leading to v\\ = 1.738 ± 0.005, in a good 
agreement with v\\ for DP in two-dimensions [Q. 

The calculation of avalanche exponents for our model 
reduces to the problem of determining the distribution 
of cluster-sizes of surface clusters in a directed site-bond 
percolation model where the concentration of bonds has 
a power-law profile (^ . In the renormalization-group ap- 
proach X = 1/z^ii is a marginal case and the cluster 
exponents may depend on the amplitude A. 

Let G(P_L,i?||) be the probability that the site 
(Px,P||) topples if a particle is added at (0,0) in the 
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steady state. Since p*{p) — p{i,p) <C 1 for large £, we can 
show that to leading order of perturbation 



G(i?i,i?,l) « Go(i?±,i?||) exp 



di/m 



(4) 



where Go(^J-i^||) is the two-point correlation function 
for the critical DP process. Using Eq. (0) we get 



G(i?±,i?||) = GaiR±,R\\)Rn 



(5) 



The value of B selected by the steady state is determined 
by the requirement that the average outflux of parti- 
cles per avalanche from the i?||-th layer equals one, i.e., 
J2r^^ G{R±,R\i) ^ 1. Since in DP the average outflux is 



J2r, Go(i?_L,i?[| 



R 



[(<i-l)yi-2/3]/z.|l 



, it follows that 
2/3]/i.|| , (6) 



where and vj^ are as above the DP exponents. 

Thus both the power-law tail and the amplitude B are 
expressed in terms of standard DP exponents. These 
facts, in turn, determine the statistics of avalanche clus- 
ters. In addition to the exponents for the distributions of 
duration, tj, and size, r^, we also define the anisotropy 
exponent ( for the average transverse extent R± ~ £^ 
of a cluster of length £. Near the DP critical line R± is 
expected to have the scaling behavior as 
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Notice that due to the power-law profile (|^), the argu- 
ment of the scaling function (f> in Eq. remains finite 
in the limit £ c», and thus C = (dp — ^±/i^\\- 

In the critical DP, the probability that a perturbation 
survives up to layer T varies as Po{T) ^ T^^/^n . As an 
expression similar to Eq. (||) applies also to the survival 
probability in the steady state of our model, we have that 
P{T) = Po{T)T-^. Using B from Eq. (i) this gives 



1 + (d - 1)C - l3/iy» 



(8) 



Using standard scaling arguments for the directed SOC 
system we notice first that the expected number of 
topplings in a cluster of length £ scales as <s>i ~ £'^*, 
that is, = Tt. Together with (t^ — = Tj — 1 this 
then gives the renormalized Ts exponent as 



= 2 - i/n 



(9) 



Inserting the best known numerical values of the expo- 
nents for two-dimensional DP gives Tt— 1.47244, Tg 
= 1.32059, and C=0.63261 . We have checked these pre- 
dictions against numerical simulations of the exponents 
and fractal dimension _D||. In the inset to Fig. 2 various 
scaling exponents are plotted versus p in the scaling re- 
gion p > p* ■ Away from a small crossover region near the 



point p = 1, the obtained values of the exponents are in- 
dependent of p within numerical error. We find Tt =1.460 
±0.014, T,= 1.313 ± 0.012, and C= 0.624 ±0.014 in fair 
agreement with the above conclusions. 

For d =3 using Eqs. (|[]^) and known numerical values 
of DP exponents ^ we find n =1.674 and =1.403. 
The upper critical dimension of our stochastic model is 
dc =5, in contrast to dc =3 in the deterministic limit 
p =1. For d >5, the DP critical exponents are /3 =1, 
=1, and iy± =1/2, leading to B =0, and thus the 
exponents have the mean-field values Tt =2, and Ts=3/2. 

In conclusion, we have demonstrated that nearness of 
the steady states to the directed-percolation critical line 
and the conservation of number of particles in the bulk 
are responsible for the emergent spatial structures in our 
stochastic sandpile model. A power-law density profile 
has been found and the self-organized criticality which 
is in a different universality class from the determinis- 
tic limit. In all dimensions d the scaling exponents of 
avalanches have been determined in terms of standard 
directed percolation critical exponents. 
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FIG. 1. Example of an avalanche running from left to 
right (dark) for p — p* on the lattice of linear size L=128 . 
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FIG. 2. Double logarithmic plot of the integrated distri- 
bution of durations P{T,p) vs. T for L=200, and p=l, 0.8, 
0.70548, 0.69, 0.68, and 0.65 (top to bottom). First two curves 
have been shifted vertically for easier display. Inset: Scaling 
exponents: (o) (O) n - 1, (A) - 1, and (★) (r^ - 1) 
plotted against p in the scaling region. 
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FIG. 3. Plot of log D{s, L) vs. log s for p = p* and for 
three different lattice sizes L=50, 100, and 200 . Inset: 
Data collapse according to Eq. (|^), where we used the val- 
ues Tt - 1=0.45 and -D||=1.45 . 
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FIG. 4. Density profile p{l) for p = p* plotted against 
the distance £ from the top (triangles) and the theoretical 



curve p{£) 



0.39 



(full line). Inset: Data in dou- 



ble-logarithmic plot. The slope is x = 0.575 ± 0.005. Within 
numerical error x remains p-independent in the region p > p* ■ 
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